arXiv: 1506.04403vl [cond-mat.mtrl-sci] 14Jun2015 


Sculpting the band gap: a computational approach 

Kiran Prasai, 1 Parthapratim Biswas, 2 and D. A. Drabold 1 

1 Department of Physics and Astronomy, Clippinger Laboratories 
Ohio University, Athens, OH 45701 
2 Department of Physics and Astronomy 
The University of Southern Mississippi, Hattiesburg, MS 39406 
(Dated: June 16, 2015) 

Materials with optimized band gap are needed in many specialized applications. In this work, 
we demonstrate that Hellmann-Feynman forces associated with the gap states can be used to find 
atomic coordinates with a desired electronic density of states. Using tight-binding models, we show 
that this approach can be used to arrive at electronically designed models of amorphous silicon 
and carbon. We provide a simple recipe to include a priori electronic information in the formation 
of computer models of materials, and prove that this information may have profound structural 
consequences. An additional example of a graphene nanoribbon is provided to demonstrate the 
applicability of this approach to engineer 2-dimensional materials. The models are validated with 
plane-wave density functional calculations. 


The central goal of materials science is the develop¬ 
ment of materials with novel properties. In general, this 
program of materials engineering has proceeded largely 
by experimental exploration. In this Letter, we offer a 
novel and direct approach to determining structures (e.g. 
atomic coordinates) that yield desired electronic or op¬ 
tical properties. The method is direct in that an ini¬ 
tial structure is purposefully modified to push the model 
toward a desired electronic density of states (for exam¬ 
ple, engineering a gap of desired magnitude, eliminating 
defect states in the gap, or perhaps changing the struc¬ 
ture of band tails, a serious issue in some photo voltaic 
(PV) applications [l]]. Such a tool is especially valuable 
for semiconductors, and may have value for chalcogenide 
phase change memory materials and device applications. 
The outcome of these simulations is a set of coordinates 
that satisfy desired optical properties, within the accu¬ 
racy of the Hamiltonian employed. Naturally, there is 
no a priori guarantee that a model with preferred elec¬ 
tronic properties should be a minimum of the total energy 
functional, but we show in two examples how to find elec¬ 
tronically optimized configurations that are indeed sta¬ 
ble, and that the method is a powerful aid to finding 
realistic models of a-Si using an electronically modified 
liquid-quench approach (normally slow and impractical, 
but effective in our approach). We provide a means to 
obtain relaxed models of a -C with widely varied sp 2 /sp 3 
ratios, and employ the method to show how to open up 
a large optical gap in a graphene nanoribbon. Models 
are validated with plane-wave density functional theory 

(dft)@,I3. 

To introduce our method, we adopt tight-binding 
Hamiltonians and employ Hellmann-Feynman forces [J, 
[|] in a novel way to determine structures with desired 
optical gap. Recall that the spatially non-local part of 
the interatomic force, the band-structure force, has the 
form: 
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Here i indexes eigenvalues (or bands), R a are the 3 N po¬ 
sitional degrees of freedom, H is the Hamiltonian, and 
is an eigenvector. For a complete basis set (such as plane 
waves), the contributions from the nuclear derivatives of 
the wave function in the first two terms cancel exactly 
(the Hellmann-Feynman theorem SI)- If one consid¬ 
ers individual terms in the sum in Eq. ©. the term FF s 
represents the contribution from the i th band (or eigen¬ 
value) to the total bandstructure force. In effect, F^ s 
is a gradient for the i th energy eigenvalue As such, 
Ff s provides the direction in the 3N-dimensional config¬ 
uration space of most rapid change of Aj. Thus, to shift 
A i to higher (lower) energies, we should move atoms in¬ 
crementally along the direction -Ff s ( + Ff s ). So, we 
simply move the atoms along (or against) these gradients 
to push tail/defect states away from the gap into the va¬ 
lence or conduction band. For small displacements 5R a 
along this gradient, the shift SAi of an eigenvalue A$ can 
be written as, 

= -3fa SR o ■ ( 2 ) 
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To this end, we introduce the term gap force for state i 
to indicate the force (negative nuclear gradient) associ¬ 
ated with eigenvalue Aj. We exploit such forces to push 
eigenvalues out of a spectral range that we wish to be 
free of states. In semiconductors or insulators, the gap 
region may be sparsely populated with defect and im¬ 
purity states, possibly localized, in which case the gap 
forces can be viewed as local perturbations arising from 
an effective energy functional Q, 
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The sum in the last term in eq. J3]) is restricted to gap 
states we wish to clear in a spectral range [£ mm , E max \. 
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Here, g(X n ) = +1 or -1 for \ n > A homo or A„ < A homo 
respectively, and /, is the occupation number of i th en¬ 
ergy level, which is either 0, 1, or 2. The parameter 7 
controls the strength of the gap force, £/ is the Fermi 
energy, and U r is the repulsive ion-ion interaction. The 
force associated with the a th degree of freedom is given 

by, 
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which can be used to obtain strong local minima by min¬ 
imizing the total energy and forces via MD simulations 
and/or relaxations. In the tight-binding formulation, the 
forces on the right-hand side of (Eq.[4]) are: 
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We show empirically that the method works well even 
for midgap states near ep. We have observed that the 
method is also applicable in the opposite mode: to maxi¬ 
mize the density of states at the Fermi level by shepherd¬ 
ing eigenvalues toward the Fermi level [7]. This might, for 
example, introduce new structural features and produce 
models with interesting electrical conductivity. 

In the rest of the paper, the first two examples pro¬ 
vide new equilibrium (i.e. relaxed) models of a-Si and 
a- C using biased dynamics to obtain a desired gap. For 
these, we perform melt-quench simulations with biased 
forces as in Eq. [I] and always take 7=1. Since we seek 
models at a minimum of the unmodified tight-binding 
total-energy functional (Eq. [3] with 7 = 0), we then re¬ 
lax the system by damping the velocity of atoms until 
the forces on atoms vanish. This simplest prescription 
has proven adequate to create with the additional flexi¬ 
bility of “optical sculpting”. The third example, tersely 
described, reports a structure of a graphene nanoribbon 
without concern for making physical forces vanish. 

We undertake our first calculation on a-Si. It is known 
for a-Si, that rapid melt-quenching from the liquid pro¬ 
duces many gap states, remnant of the 6-fold liquid metal 
and an unsatisfactory model for a tetrahedral amorphous 
semiconductor. In this example, we impose biased dy¬ 
namics favoring the creation of a gap in the range ob¬ 
served in the best available WWW models Q. After equi¬ 
librating the liquid in the conventional way, we quench 
with dissipative dynamics (velocity rescaling to 300 K) 
and biased forces. We note, for the choice of 7 = 1 we 
invoke, that the average gap forces remain less than 20% 
of the TBMD forces on atoms: the dynamics is modi¬ 
fied in a somewhat subtle way, but operating over many 
steps, the method yields structures improved both opti¬ 
cally (by construction) but also structurally , a consider¬ 
able bonus and proof that inclusion of electronic a pri¬ 
ori information does influence structure, and in a way 


that improves agreement with experiments. We use the 
Goodwin-Skinner-Pettifor (GSP) Hamiltonian for Si @. 

In the quenching process, we add gap forces to push the 
eigenvalues away from the Fermi level, using prescription 
U Since we relax the model with 7 = 0 after quenching, 
the final model is at equilibrium according to the true 
forces derived from [3] with 7 = 0. 

A 216-atom model of liquid silicon (l- Si) was prepared 
by cooling an initial random configuration from a temper¬ 
ature of 2500 K to the melting point « 1780 K in several 
steps, which were followed by equilibration for a period 
of 50 ps per step and total-energy relaxations. We have 
verified that our l -Si model produces features similar to l- 
Si models obtained from the GSP Hamiltonian by earlier 
workers [lOj-QjJ. In figs. [HU we present structural and 
electronic properties from the dynamics for 7 = 0 (e.g. 
conventional MD) and 1 starting with liquid silicon. 

Figure [T] compares the radial and bond-angle distri¬ 
bution functions for three models of a-Si obtained from 
conventional TBMD, biased TBMD, and the WWW 
method Q, the last being the “gold standard” of a-Si 
models. A comparison of the RDFs and an analysis 
of bond lengths shows that the biased dynamics signifi¬ 
cantly increases the short-ranged order: the biased model 
exhibits fewer short and long bonds as compared to con¬ 
ventional TBMD and WWW models. A remarkable fea¬ 
ture of the biased-MD model is the absence of 60° bond 
angles. These angles are typically associated with frozen 
liquid-like configurations, which plague conventional MD 
simulations by producing 3-fold defects (dangling bonds). 
The latter are notoriously difficult to remove from MD 
models via conventional means and have detrimental ef¬ 
fects on electronic and optical properties. It is notewor¬ 
thy that the inclusion of gap forces can eliminate these 
unphysical features completely and impart more tetra¬ 
hedral order in the structure. Atomic coordination is 
also markedly improved in biased dynamics. Atomic- 
coordination statistics show that 97.2% atoms in the 
biased-MD model are 4-fold coordinated. This result is 
not only superior to 87% 4-fold coordination in conven¬ 
tional MD but also better than earlier works rep orted in 
the literature using the GSP Hamiltonian [Io[ [lj, [l3| . 
The biased-TBMD model has fewer defects around the 
Fermi-energy than the TBMD model (fig. [2]). Next, 
we check these models using using DFT in local den¬ 
sity approximation (LDA). LDA calculations with VASP 
show that the energy of the biased-TBMD configuration 
is lower than that of the conventional TBMD configu¬ 
rations. Imposing electronic constraints leads to relaxed 
models in better agreement with structural experiments. 

To explore the reproducibility of our method, we sam¬ 
pled 25 well-separated l -Si configurations as starting co¬ 
ordinates and quenched these configurations with 7 = 
0 and 1 as described earlier. This was followed by total- 
energy relaxations of the models to their respective local 
minimum. We analyzed these models to obtain the num¬ 
ber of atoms with 4-fold coordination as a measure of the 
merit of the model. Of 25 biased-TBMD models that we 






FIG. 1: (Color online) Radial (left) and bond-angle (right) 
distributions of a-Si from three models: biased TBMD (red), 
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FIG. 2: (Color online) Density of electronic states (EDOS) 
near the band-gap region for three models of a-Si: TBMD, 
biased TBMD and WWWQ. The density of states within 
the gap is shown as a vertical line with respective color as 
indicated. Inset shows full EDOS. 


studied in this work, 20 models were found to produce 
better (more tetrahedral) topology than TBMD models. 

As a second example we study a-C. Tetrahedral amor¬ 
phous carbon (ta-C), has some properties reminiscent of 
diamond while potentially holding some advantages [Til 
im The tight-binding model of Xu et al [T(| has been 
used previously to model ta-C with limited success G3- 
These calculations involved a quench from a high-density 
liquid (/-C) and volume rescaling at lower temperature. 
Using the same Hamiltonian, we demonstrate that a sim¬ 
pler melt-quench method can yield improved models. 
Amorphous carbon dominated by sp 3 bonding is char¬ 
acterized by a band gap of about 2 eV (depending on the 
fraction of sp 3 bonded atoms) in contrast to sp 2 -bonded 
a-C, which has a gap of less than 0.5 eV |l4|. The per¬ 
fectly sp 3 -bonded WWW model of ta-C j§j, relaxed with 
the Xu Hamiltonian, has a gap of 4.1 eV. We used this 
spectral range with biased melt-quench TBMD to form 
models without states in the gap region as exhibited by 
WWW models. 

Starting with liquid carbon of density 3.5 gm/cm 3 
equilibrated at 10000 K, we quenched the model to 700 



FIG. 3: (Color online) Radial (left) and bond-angle (right) 
distribution functions of ta-C from biased TBMD, TBMD, 
and WWW. See text for discussion of results. 


K at a rate of 500 K/ps. At this point, two parallel runs 
(quenching) were performed: one with regular TBMD 
forces (7 — 0) and the other with additional gap forces 
for 7=1. Following first example, we designate these 
runs as ‘TBMD’ and ‘biased TBMD’, respectively. Both 
quenched models were relaxed to their respective local 
minimum until the force on each atom is less than 0.05 
eV/A. Gap forces were switched off during post-quench 
relaxation for biased-TBMD models. The magnitude of 
gap forces were found to be smaller (< 20%) than the 
corresponding TBMD forces throughout the quenching 
process. 

The unaided TBMD with the Xu Hamiltonian prefers 
sp 2 dominated network as observed in our calculations 
and in Refs. ia The diamond-like sp 3 -bonded net¬ 
works reported in m appear to be an artifact of high 
density or high pressure on l- C. Our calculations pro¬ 
duced models with up to 94% 4-fold coordination com¬ 
pared with 74% and 89% in E3- Our results are readily 
reproducible and do not involve arbitrary manipulation 
of density. We have conducted 25 quenching runs with 
different starting liquid models and all of these models 
produce tetrahedral networks with more than 90% 4-fold 
coordination. The structural features, including RDF, 
BDF and atomic coordinations from biased TBMD re¬ 
semble closely with ta-C WWW model. The TBMD 
model is dominated by sp 2 bonding and registers distinct 
peaks in the RDF and BDF (Figure [3]). The density of 
electronic states of the biased-TBMD model shows that 
the gap opens up to 0.7 eV as compared to 0.21 eV in the 
TBMD model. Also, the biased-TBMD model has only 
14 states in the gap region exhibited by WWW mod¬ 
els as compared to 71 states in the TBMD model: see 
Fig. [H The electronic structure of these models is also 
confirmed by LDA calculations using VASP. The biased- 
TBMD model shows few scattered states in the gap re¬ 
gion as opposed to the ‘metal-like’ electronic structure of 
the TBMD model. Total energy calculations using LDA 
show that the biased-TBMD model has 0.31 eV/atom 
less energy than the regular TBMD model. The biased 
model is also stable under relaxation using LDA. Such 
relaxation decreases the total energy by 0.07 eV/atom 
while preserving structural ordering of the model. 
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FIG. 4: (Color online) Density of electronic states (EDOS) 
for three models of a-C: TBMD, biased TBMD and 
WWW[8,]. The upper plot is obtained using the 
tight-binding model of [lj| and the lower one from LDA 

an 


Finally, we have applied this approach to graphene 
nanoribbons. Using a T-point calculation with the same 
TB Hamiltonian in example 2 on an armchair graphene 
nano-ribbon (AGNR) of width 11.07 A, we have carried 
out low-temperature MD simulations biased toward an 
enlarged bandgap, without trying to drive the model to 
a minimum of the total energy. Depending on the value 
of 7 selected, we obtain structures with gaps of up to 1.58 
eV (wider gap is observed using LDA). But these com¬ 
puter models have high strain on the edge atoms and 
the length of C-C dimers along the edge becomes smaller 
than the average bond length by 18%. 

As we have shown, our method is best employed in 
a “statistical mode”-unsurprisingly the final structures 
depend on the initial state. In some fraction (s=s 20%) 
the method does not improve the gap in case of a-Si. We 
suppose that this may be due to the very simple rule of 


shifting atoms along gradients toward the nearest band 
edge, even for eigenvalues very near Ef. 

For a-Si, we use a priori knowledge of the gap from the 
best available models. In the general case, one can define 
a gap by trial and error, with the choice being determined 
in part by a requirement that the physical forces vanish 
at the end. For a-C, considerable flexibility is afforded 
by our approach in tuning sp 1 2 / sp 3 4 5 6 7 * 9 ratios. We have ex¬ 
perimented with various 7 , and have found no particular 
advantage to selecting 7 ^ 1 to date. Furthermore, pre¬ 
liminary studies suggest that the results presented here 
also accrue for larger (512-atom) models. We expect that 
the scheme will be useful for many other complex mate¬ 
rials not only for discovering structures with desired gaps 
but also for imposing electronic constraints in modeling. 

As is the case with all methods, our approach has lim¬ 
itations: (1) For this first report we use standard tight- 
binding Hamiltonians for the simulations. Such Hamil¬ 
tonians are well known to have imperfect transferability 
(for this reason we are currently extending the scheme to 
plane-wave DFT, a straightforward but tedious under¬ 
taking) and 2 ) even in a density-functional framework, 
gap estimates from Kohn-Sham eigenvalues are spurious, 
though usually these account reasonably well for trends. 
With significant computational expense, these estimates 
may be improved e.g. with GW or Hybrid Functional 
schemes [ij. Despite these limitations, we demonstrate 
the utility of the method with two examples and suggest 
that the approach may be developed in promising ways. 
Nevertheless, it is plain that the method is an extremely 
useful new tool, even using this simplest implementation. 
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